library(rgrass7)
## Loading required package: XML
## GRASS GIS interface loaded with GRASS version: GRASS 7.8.2 (2019)
## and location: pantuflas
gisdbase <- 'grass-data-test' #Base de datos de GRASS GIS
wd <- getwd() #Directorio de trabajo
wd
## [1] "/home/jr/unidad-4-asignacion-1-procesos-fluviales"
loc <- initGRASS(gisBase = "/usr/lib/grass78/",
home = wd,
gisDbase = paste(wd, gisdbase, sep = '/'),
location = 'pantuflas',
mapset = "PERMANENT",
override = TRUE)
gmeta()
## gisdbase /home/jr/unidad-4-asignacion-1-procesos-fluviales/grass-data-test
## location pantuflas
## mapset PERMANENT
## rows 892
## columns 971
## north 2122920
## south 2042640
## west 253800
## east 341190
## nsres 90
## ewres 90
## projection +proj=utm +no_defs +zone=19 +a=6378137 +rf=298.257223563
## +towgs84=0.000,0.000,0.000 +type=crs +to_meter=1
execGRASS(
'g.list',
flags = 't',
parameters = list(
type = c('raster', 'vector')
)
)
## raster/demint
## raster/rbasin_pant_demint_accumulation
## raster/rbasin_pant_demint_aspect
## raster/rbasin_pant_demint_dist2out
## raster/rbasin_pant_demint_drainage
## raster/rbasin_pant_demint_hack
## raster/rbasin_pant_demint_hillslope_distance
## raster/rbasin_pant_demint_horton
## raster/rbasin_pant_demint_shreve
## raster/rbasin_pant_demint_slope
## raster/rbasin_pant_demint_strahler
## raster/stream-de-rstr
## vector/rbasin_pant_demint_basin
## vector/rbasin_pant_demint_mainchannel
## vector/rbasin_pant_demint_network
## vector/rbasin_pant_demint_outlet
## vector/rbasin_pant_demint_outlet_snap
## vector/stream_de_rstr
Convertir a números enteros la extensión y la resolución del DEM
library(raster)
## Loading required package: sp
rutadem <- 'data/dem.tif'
rawextent <- extent(raster(rutadem))
rawextent
## class : Extent
## xmin : 253740.1
## xmax : 341262.4
## ymin : 2042609
## ymax : 2122949
devtools::source_url('https://raw.githubusercontent.com/geofis/rgrass/master/integerextent.R')
## SHA-1 hash of file is 65b57fe7cb20cd1976572cd0a7e98def9e95d83c
devtools::source_url('https://raw.githubusercontent.com/geofis/rgrass/master/xyvector.R')
## SHA-1 hash of file is 89fa5ae436d9a7d7a0c799b789c560eb5e421cfd
newextent <- intext(e = rawextent, r = 90, type = 'inner')
newextent
## class : Extent
## xmin : 253800
## xmax : 341190
## ymin : 2042640
## ymax : 2122920
gdalUtils::gdalwarp(
srcfile = 'data/dem.tif',
dstfile = 'data/demint.tif',
te = xyvector(newextent),
tr = c(90,90),
r = 'bilinear',
overwrite = T
)
## Registered S3 method overwritten by 'R.oo':
## method from
## throw.default R.methodsS3
## NULL
Importar a sesión de GRASS
rutademint <- 'data/demint.tif'
execGRASS(
"g.proj",
flags = c('t','c'),
georef=rutademint)
gmeta()
execGRASS(
"r.in.gdal",
flags='overwrite',
parameters=list(
input=rutademint,
output="demint"
)
)
execGRASS(
"g.region",
parameters=list(
raster = "demint",
align = "demint"
)
)
gmeta()
## gisdbase /home/jr/unidad-4-asignacion-1-procesos-fluviales/grass-data-test
## location pantuflas
## mapset PERMANENT
## rows 892
## columns 971
## north 2122920
## south 2042640
## west 253800
## east 341190
## nsres 90
## ewres 90
## projection +proj=utm +no_defs +zone=19 +a=6378137 +rf=298.257223563
## +towgs84=0.000,0.000,0.000 +type=crs +to_meter=1
execGRASS(
'g.list',
flags = 't',
parameters = list(
type = c('raster', 'vector')
)
)
## raster/demint
## raster/rbasin_pant_demint_accumulation
## raster/rbasin_pant_demint_aspect
## raster/rbasin_pant_demint_dist2out
## raster/rbasin_pant_demint_drainage
## raster/rbasin_pant_demint_hack
## raster/rbasin_pant_demint_hillslope_distance
## raster/rbasin_pant_demint_horton
## raster/rbasin_pant_demint_shreve
## raster/rbasin_pant_demint_slope
## raster/rbasin_pant_demint_strahler
## raster/stream-de-rstr
## vector/rbasin_pant_demint_basin
## vector/rbasin_pant_demint_mainchannel
## vector/rbasin_pant_demint_network
## vector/rbasin_pant_demint_outlet
## vector/rbasin_pant_demint_outlet_snap
## vector/stream_de_rstr
Generar red de drenaje para obtener coordenada posteriormente
execGRASS(
"r.stream.extract",
flags = c('overwrite','quiet'),
parameters = list(
elevation = 'demint',
threshold = 80,
stream_raster = 'stream-de-rstr',
stream_vector = 'stream_de_rstr'
)
)
## Warning in execGRASS("r.stream.extract", flags = c("overwrite", "quiet"), : The command:
## r.stream.extract --overwrite --quiet elevation=demint threshold=80 stream_raster=stream-de-rstr stream_vector=stream_de_rstr
## produced at least one warning during execution:
## WARNING: Vector map <stream_de_rstr> already exists and will be overwritten
## WARNING: Vector map <stream_de_rstr> already exists and will be overwritten
execGRASS(
'g.list',
flags = 't',
parameters = list(
type = c('raster', 'vector')
)
)
## raster/demint
## raster/rbasin_pant_demint_accumulation
## raster/rbasin_pant_demint_aspect
## raster/rbasin_pant_demint_dist2out
## raster/rbasin_pant_demint_drainage
## raster/rbasin_pant_demint_hack
## raster/rbasin_pant_demint_hillslope_distance
## raster/rbasin_pant_demint_horton
## raster/rbasin_pant_demint_shreve
## raster/rbasin_pant_demint_slope
## raster/rbasin_pant_demint_strahler
## raster/stream-de-rstr
## vector/rbasin_pant_demint_basin
## vector/rbasin_pant_demint_mainchannel
## vector/rbasin_pant_demint_network
## vector/rbasin_pant_demint_outlet
## vector/rbasin_pant_demint_outlet_snap
## vector/stream_de_rstr
Obtener coordenada
library(sp)
use_sp()
library(mapview)
netw <- spTransform(
readVECT('stream_de_rstr'),
CRSobj = CRS("+init=epsg:4326"))
mapview(netw, col.regions = 'blue', legend = FALSE)
## Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
## deprecated. It might return a CRS with a non-EPSG compliant axis order.
## Warning: sf layer has inconsistent datum (+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs).
## Need '+proj=longlat +datum=WGS84'
Ejecutar r.basin
pref <- 'rbasin_pant'
execGRASS(
"r.basin",
flags = 'overwrite',
parameters = list(
map = 'demint',
prefix = pref,
coordinates = outlet,
threshold = 80,
dir = 'salidas-rbasin/pantuflas'
)
)
execGRASS(
'g.list',
flags = 't',
parameters = list(
type = c('raster', 'vector')
)
)
## raster/demint
## raster/rbasin_pant_demint_accumulation
## raster/rbasin_pant_demint_aspect
## raster/rbasin_pant_demint_dist2out
## raster/rbasin_pant_demint_drainage
## raster/rbasin_pant_demint_hack
## raster/rbasin_pant_demint_hillslope_distance
## raster/rbasin_pant_demint_horton
## raster/rbasin_pant_demint_shreve
## raster/rbasin_pant_demint_slope
## raster/rbasin_pant_demint_strahler
## raster/stream-de-rstr
## vector/rbasin_pant_demint_basin
## vector/rbasin_pant_demint_mainchannel
## vector/rbasin_pant_demint_network
## vector/rbasin_pant_demint_outlet
## vector/rbasin_pant_demint_outlet_snap
## vector/stream_de_rstr
Si r.basin arrojara error (sólo en el caso de error, no en caso de advertencia), ejecutar este bloque para borrar las salidas anteriores y reejecutar el r.basin:
execGRASS(
"g.remove",
flags = 'f',
parameters = list(
type = c('raster','vector'),
pattern = paste0(pref, '*')
)
)
Explorar los parámetros de cuenca
library(readr)
rbpantpar1 <- read_csv("salidas-rbasin/pantuflas/rbasin_pant_demint_parametersT.csv")
## Parsed with column specification:
## cols(
## .default = col_double(),
## Rectangle_containing_basin_N_W = col_character(),
## Rectangle_containing_basin_S_E = col_character()
## )
## See spec(...) for full column specifications.
rbpantpar1 %>% tibble::as_tibble()
## # A tibble: 1 x 34
## x y Easting_Centroi… Northing_Centro… Rectangle_conta…
## <dbl> <dbl> <dbl> <dbl> <chr>
## 1 313155 2.09e6 319815 2087685 ('313110', '209…
## # … with 29 more variables: Rectangle_containing_basin_S_E <chr>,
## # Area_of_basin_km2 <dbl>, Perimeter_of_basin_km <dbl>,
## # Max_Elevation <dbl>, Min_Elevation <dbl>, Elevation_Difference <dbl>,
## # Mean_Elevation <dbl>, Mean_Slope <dbl>,
## # Length_of_Directing_Vector_km <dbl>,
## # Prevalent_Orientation_deg_from_north_ccw <dbl>,
## # Compactness_Coefficient <dbl>, Circularity_Ratio <dbl>,
## # Topological_Diameter <dbl>, Elongation_Ratio <dbl>,
## # Shape_Factor <dbl>, Concentration_Time_hr <dbl>,
## # Length_of_Mainchannel_km <dbl>,
## # Mean_slope_of_mainchannel_percent <dbl>,
## # Mean_hillslope_length_m <dbl>, Magnitudo <dbl>,
## # Max_order_Strahler <dbl>, Number_of_streams <dbl>,
## # Total_Stream_Length_km <dbl>, First_order_stream_frequency <dbl>,
## # Drainage_Density_km_over_km2 <dbl>, Bifurcation_Ratio_Horton <dbl>,
## # Length_Ratio_Horton <dbl>, Area_ratio_Horton <dbl>,
## # Slope_ratio_Horton <dbl>
rbpantpar2 <- read_csv(
"salidas-rbasin/pantuflas/rbasin_pant_demint_parameters.csv",
skip=2, col_names = c('Parameter', 'Value'))
## Parsed with column specification:
## cols(
## Parameter = col_character(),
## Value = col_character()
## )
rbpantpar2 %>% print(n=Inf)
## # A tibble: 32 x 2
## Parameter Value
## <chr> <chr>
## 1 Easting Centroid of basin 319815.00
## 2 Northing Centroid of basin 2087685.00
## 3 Rectangle containing basin N-W ('313110', '20975…
## 4 Rectangle containing basin S-E ('326880', '20782…
## 5 Area of basin [km^2] 128.6087625
## 6 Perimeter of basin [km] 63.9170486606335
## 7 Max Elevation [m s.l.m.] 2619.45
## 8 Min Elevation [m s.l.m.] 1054.484
## 9 Elevation Difference [m] 1564.966
## 10 Mean Elevation 1588.877
## 11 Mean Slope 13.74
## 12 Length of Directing Vector [km] 7.367794853278693
## 13 Prevalent Orientation [degree from north, counterclo… 0.441915859817346…
## 14 Compactness Coefficient 4.994895129597467
## 15 Circularity Ratio 0.395591541103538…
## 16 Topological Diameter 18.0
## 17 Elongation Ratio 0.6185502219049736
## 18 Shape Factor 6.216632397877906
## 19 Concentration Time (Giandotti, 1934) [hr] 2.413889375890686
## 20 Length of Mainchannel [km] 20.687850635
## 21 Mean slope of mainchannel [percent] 6.914012642018835
## 22 Mean hillslope length [m] 858.9935
## 23 Magnitudo 33.0
## 24 Max order (Strahler) 4
## 25 Number of streams 44
## 26 Total Stream Length [km] 99.8046
## 27 First order stream frequency 0.2565921587185787
## 28 Drainage Density [km/km^2] 0.7760326595164928
## 29 Bifurcation Ratio (Horton) 3.2429
## 30 Length Ratio (Horton) 0.8840
## 31 Area ratio (Horton) 3.7251
## 32 Slope ratio (Horton) 1.5582
Limpiar archivo de bloqueo del conjunto de mapas de GRASS
unlink_.gislock()
Referencias